# ============================================================================
# 线弹性GPU测试 - 验证GPU求解器基础功能
# 先验证线弹性，再考虑非线性
# ============================================================================

wipe
puts "线弹性GPU vs CPU对比测试"

# ============================================================================
# 模型配置
# ============================================================================

model basic -ndm 2 -ndf 2

# 优化网格设计
set nEleX 50
set nEleY 50
set dx 0.1
set dy 0.1

puts "创建线弹性测试模型: ${nEleX}x${nEleY} = [expr $nEleX*$nEleY] 个单元"

# 创建节点
set nodeTag 1
for {set j 0} {$j <= $nEleY} {incr j} {
    for {set i 0} {$i <= $nEleX} {incr i} {
        node $nodeTag [expr $i * $dx] [expr $j * $dy]
        incr nodeTag
    }
}

# ============================================================================
# 选择监测节点用于CPU vs GPU结果对比
# ============================================================================

# 选择顶部中心节点作为监测点
set monitorNode [expr $nEleY*($nEleX+1) + ($nEleX+1)/2]
puts "监测节点: $monitorNode (顶部中心节点)"

# ============================================================================
# 线弹性材料模型 - 从DruckerPrager参数转换
# ============================================================================

# 从 K=100MPa, G=75MPa 计算 E 和 nu
# E = 9KG/(3K+G) = 9*100*75/(300+75) = 180 MPa
# nu = (3K-2G)/(2(3K+G)) = (300-150)/(2*375) = 0.2

set E 180.0e6    ;# 杨氏模量 180 MPa
set nu 0.2       ;# 泊松比
set rho 1800.0   ;# 密度

nDMaterial ElasticIsotropic 1 $E $nu $rho

puts "材料: 线弹性各向同性 E=$E Pa, nu=$nu"

# 创建线弹性四边形单元
set eleTag 1
for {set j 0} {$j < $nEleY} {incr j} {
    for {set i 0} {$i < $nEleX} {incr i} {
        set n1 [expr $j*($nEleX+1) + $i + 1]
        set n2 [expr $j*($nEleX+1) + $i + 2]
        set n3 [expr ($j+1)*($nEleX+1) + $i + 2]
        set n4 [expr ($j+1)*($nEleX+1) + $i + 1]

        # 使用平面应力
        element quad $eleTag $n1 $n2 $n3 $n4 1.0 "PlaneStress" 1
        incr eleTag
    }
}

# ============================================================================
# 边界条件和加载
# ============================================================================

# 底部固定边界
for {set i 1} {$i <= [expr $nEleX+1]} {incr i} {
    fix $i 1 1
}

# 侧边边界：允许竖向位移，限制水平位移
for {set j 1} {$j <= $nEleY} {incr j} {
    set leftNode [expr $j*($nEleX+1) + 1]
    set rightNode [expr $j*($nEleX+1) + ($nEleX+1)]
    fix $leftNode 1 0
    fix $rightNode 1 0
}

# 线性加载 - 适度荷载
pattern Plain 1 Linear {
    # 顶部施加分布荷载 1000 N/m² (线弹性不需要太大)
    for {set i 1} {$i <= [expr $nEleX+1]} {incr i} {
        set topNode [expr $nEleY*($nEleX+1) + $i]
        load $topNode 0.0 -1000.0
    }
}

set totalDOF [expr ($nEleX+1)*($nEleY+1)*2]
puts "模型: [expr $nEleX*$nEleY] 单元, ${totalDOF} 自由度"
puts "预期：线弹性问题应该1步收敛"

# ============================================================================
# CPU 线弹性分析
# ============================================================================

puts "\n=== CPU 线弹性分析 ==="

constraints Plain
numberer RCM
system BandGeneral
test NormDispIncr 1.0e-8 10  ;# 线性问题应该很快收敛
algorithm Newton
integrator LoadControl 1.0  ;# 一步加载
analysis Static

# 删除已存在的输出文件
file delete -force "cpu_response.txt"

# 创建记录器
recorder Node -file "cpu_response.txt" -time -node $monitorNode -dof 1 2 disp

puts "开始CPU分析..."
set cpuStart [clock clicks -milliseconds]
set ok [analyze 1]
set cpuEnd [clock clicks -milliseconds]
set cpuTime [expr $cpuEnd - $cpuStart]

if {$ok == 0} {
    puts "CPU分析成功: ${cpuTime} ms"
} else {
    puts "❌ CPU分析失败!"
}

# 完全清理并重建模型
wipe

# ============================================================================
# 重建模型用于GPU分析
# ============================================================================

puts "\n重建模型用于GPU分析..."
model basic -ndm 2 -ndf 2

# 重新创建节点
set nodeTag 1
for {set j 0} {$j <= $nEleY} {incr j} {
    for {set i 0} {$i <= $nEleX} {incr i} {
        node $nodeTag [expr $i * $dx] [expr $j * $dy]
        incr nodeTag
    }
}

# 重新定义材料
nDMaterial ElasticIsotropic 1 $E $nu $rho

# 重新创建单元
set eleTag 1
for {set j 0} {$j < $nEleY} {incr j} {
    for {set i 0} {$i < $nEleX} {incr i} {
        set n1 [expr $j*($nEleX+1) + $i + 1]
        set n2 [expr $j*($nEleX+1) + $i + 2]
        set n3 [expr ($j+1)*($nEleX+1) + $i + 2]
        set n4 [expr ($j+1)*($nEleX+1) + $i + 1]
        element quad $eleTag $n1 $n2 $n3 $n4 1.0 "PlaneStress" 1
        incr eleTag
    }
}

# 重新施加边界条件
for {set i 1} {$i <= [expr $nEleX+1]} {incr i} {
    fix $i 1 1
}
for {set j 1} {$j <= $nEleY} {incr j} {
    set leftNode [expr $j*($nEleX+1) + 1]
    set rightNode [expr $j*($nEleX+1) + ($nEleX+1)]
    fix $leftNode 1 0
    fix $rightNode 1 0
}

# 重新定义荷载
pattern Plain 1 Linear {
    for {set i 1} {$i <= [expr $nEleX+1]} {incr i} {
        set topNode [expr $nEleY*($nEleX+1) + $i]
        load $topNode 0.0 -1000.0
    }
}

puts "模型重建完成"

# ============================================================================
# GPU 线弹性分析
# ============================================================================

puts "\n=== GPU 线弹性分析 ==="

constraints Plain
numberer RCM
system GPUSolver  ;# 使用GPU求解器
test NormDispIncr 1.0e-8 10
algorithm Newton
integrator LoadControl 1.0
analysis Static

# 删除已存在的输出文件
file delete -force "gpu_response.txt"

# 创建记录器
recorder Node -file "gpu_response.txt" -time -node $monitorNode -dof 1 2 disp

puts "开始GPU分析..."
set gpuStart [clock clicks -milliseconds]
set ok [analyze 1]
set gpuEnd [clock clicks -milliseconds]
set gpuTime [expr $gpuEnd - $gpuStart]

if {$ok == 0} {
    puts "GPU分析成功: ${gpuTime} ms"
} else {
    puts "❌ GPU分析失败!"
}

# ============================================================================
# 性能对比
# ============================================================================

puts "\n=========================================="
puts "线弹性GPU vs CPU性能对比"
puts "=========================================="
puts "测试模型: [expr $nEleX*$nEleY] 单元, ${totalDOF} 自由度"
puts ""

puts "CPU 性能: ${cpuTime} ms"
puts "GPU 性能: ${gpuTime} ms"
puts ""

if {$cpuTime > 0 && $gpuTime > 0} {
    set speedup [expr double($cpuTime) / $gpuTime]
    if {$speedup > 1.0} {
        puts "✓ GPU加速: [format %.2fx $speedup]"
    } else {
        puts "⚠ GPU较慢: [format %.2fx [expr 1.0/$speedup]]"
    }
}

# ============================================================================
# 结果精度对比
# ============================================================================

puts ""
puts "=========================================="
puts "CPU vs GPU 结果精度对比"
puts "=========================================="

if {[file exists "cpu_response.txt"] && [file exists "gpu_response.txt"]} {
    # 读取CPU结果
    set cpu_file [open "cpu_response.txt" r]
    set cpu_data [split [string trim [read $cpu_file]] "\n"]
    close $cpu_file

    # 读取GPU结果
    set gpu_file [open "gpu_response.txt" r]
    set gpu_data [split [string trim [read $gpu_file]] "\n"]
    close $gpu_file

    if {[llength $cpu_data] > 0 && [llength $gpu_data] > 0} {
        set cpu_line [split [lindex $cpu_data 0]]
        set gpu_line [split [lindex $gpu_data 0]]

        if {[llength $cpu_line] >= 3 && [llength $gpu_line] >= 3} {
            set x_cpu [lindex $cpu_line 1]
            set x_gpu [lindex $gpu_line 1]
            set y_cpu [lindex $cpu_line 2]
            set y_gpu [lindex $gpu_line 2]

            set diff_x [expr abs($x_cpu - $x_gpu)]
            set diff_y [expr abs($y_cpu - $y_gpu)]

            puts "监测节点位移:"
            puts [format "  X方向: CPU=%12.6e, GPU=%12.6e, 差值=%.2e" $x_cpu $x_gpu $diff_x]
            puts [format "  Y方向: CPU=%12.6e, GPU=%12.6e, 差值=%.2e" $y_cpu $y_gpu $diff_y]
            puts ""

            set max_diff [expr max($diff_x, $diff_y)]
            if {$max_diff < 1e-12} {
                puts "✓ 精度: 优秀 (差值 < 1e-12)"
            } elseif {$max_diff < 1e-10} {
                puts "✓ 精度: 很好 (差值 < 1e-10)"
            } elseif {$max_diff < 1e-8} {
                puts "✓ 精度: 良好 (差值 < 1e-8)"
            } elseif {$max_diff < 1e-6} {
                puts "⚠ 精度: 一般 (差值 < 1e-6)"
            } else {
                puts "❌ 精度: 较差 (差值 >= 1e-6)"
            }
        }
    }
} else {
    puts "❌ 找不到结果文件"
}

puts ""
puts "=========================================="
puts "测试完成!"
puts "=========================================="
puts "下一步："
puts "1. 如果精度不好，检查GPU残差计算和刚度矩阵组装"
puts "2. 如果精度很好，可以尝试加入非线性材料"
puts "=========================================="
